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Abstract 

We present a new approach to determine numerically the statistical behavior of small-scale 
structures in hydrodynamic turbulence. Starting from the functional integral representation of the 
random-force-driven Burgers equation we show that Monte Carlo simulations allow us to determine 
the anomalous scaling of high-order moments of velocity differences. Given the general applicability 
of Monte Carlo methods, this opens up the possibility to address also other systems relevant to 
turbulence within this framework. 
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The small-scale statistical properties of hydrodynamic turbulence is an old and tantalizing 
problem For turbulent flow stirred at large scales and far from the boundaries one 
expects a universal scaling for the small-scale fluctuations. Indeed, experiment gives strong 
indications for such universal behavior in Navier-Stokes turbulence jj^. The exact values 
of the scaling exponents however are still under debate. In such a situation it is useful to 
have a model system at hand that shares some essential properties with the original problem 
and allows for a clear physical understanding. 

The random-force-driven Burgers equation 

dtu + ud^u - vdlu = /(x, t) , (1) 

is one such example. It was originally conceived as a one-dimensional model for compressible 
hydrodynamic turbulence [3| and provides a useful benchmark setting to test new analytical 
and numerical methods for real-world turbulence Here, u is the velocity, and / 

a centered random field displaying Gaussian statistics. We will consider the special case 
where the system is driven by a self-similar forcing that is white in time. The two-point 
correlation function of the stochastic forcing in Fourier space is given by 

{f{k, t)f{k', t')) oc Do\kfS{k + k')S{t - t') , (2) 

where the parameter /3 determines the relative importance of the stirring mechanism at 
different scales, and the dimensionful constant Dq measures its strength. For /3 large and 
negative the forcing effectively acts at large scales. On the other hand, the kinematic 
viscosity in ([1]) provides a dissipation scale t] and for — 0"'' the two characteristic 
scales ?7, and the system size L separate. The stochastic forcing drives the system into a 
non-equilibrium steady state, where in the range 77 ^ |A;|~^ ^ L the energy flux through 



wavenumber k behaves as n(A;) oc [lOl. 



The case /3 = — 1 corresponds to the physically interesting situation where the flux n(/c) 
is constant (up to logarithmic corrections) and the interplay of the stochastic forcing and 
advective term leads to a Kolmogorov energy spectrum E{k) oc Q,!!^. The physical 

picture behind this scenario is the appearance of shocks with a finite dissipative width (see 
e.g. Fig. [1]). The large fluctuations associated with the negative gradient of the front give the 
dominant contribution to the anomalous scaling of velocity differences Am = u{x + r) — u{x). 
In particular, we have the structure functions (|Au|") oc r^", and the scaling exponents 




FIG. 1: Typical velocity profile u{x) from a simulation on a 254 x 1024 (space x time) lattice, 
where x is taken in units of the spatial lattice size L. 

Cn = 1 for n > 3 strongly deviate from the Kolmogorov scaling prediction = n/S that 
follows from a naive dimensional analysis jsl, oj. These rare fluctuations are strongly non- 
Gaussian and lead to the known asymptotic left tail of the probability distribution function 
(PDF) for velocity differences V{Au,r) jl^ . 

Here, we approach the problem from the functional integral point of view 15Nl7l|. The 



functional integral gives a non-perturbative definition of the field theory and thus, it is 
ideally suited to study the strong and rare fluctuations present in fully developed turbu- 
lence that give the main contribution to the high-order moments of velocity differences. By 
sampling the associated probability distribution functional via Monte Carlo methods we 
show that it is possible to determine the scaling behavior of structure functions from first 
principles. Monte Carlo simulations are directly transferable to other systems of interest 
and are free of any modeling assumptions. Though not directly competitive with conven- 
tional time- advancing methods as, e.g. pseudo-spectral or finite-difference methods, Monte 
Carlo simulations may provide a unique perspective on such important problems as, e.g. 



intermittency in fully developed turbulence 



scaling behavior of Burgers turbu 
mechanisms for intermittency Q, 



ence 



m 



14j |. In view of the well-established anomalous 



18| and the physical picture of the underlying 



20|, this provides an ideal setting to test our method 
and understand possible systematic effects at finite Reynolds number and system size. We 
emphasize again that in this paper we are not aiming to complete, or even improve on the 
accuracy obtained with other methods for Burgers turbulence. We rather want to provide a 
test of the generally applicable functional integral method and to demonstrate that a very 



3 



reasonable accuracy can be obtained from this approach, a fact that was highly unclear at 
the beginning of this project. 



The functional integral for t 
Martin-Siggia-Rose formalism 
We have the field theory 



le random-force-driven Burgers equation is obtained via the 
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22| by introducing an auxiliary response field fi. 
[du] [dfi]exp{—S[u, fi]} , (3) 



with the action 



+ - dtdxdyfi{x,t)D{x-y)fi{y,t) 



(4) 

where D{x — y) is the spatial part of the two-point correlation function In this form 
the action does not satisfy positivity. To obtain a Gibbs measure that can be sampled by 
a Markov chain Monte Carlo (MCMC) algorithm we integrate out the auxiliary field. This 
leaves us with the probability density functional 
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(5) 



P[u] = exp I ~ 2 J dtdxdy {dt'u + ud^u — vdlu) 
D~^{x — y){dtu + ud^u — vdlu)^ . 



The theory is then defined by placing the field u{x,t) on the sites of a regular space-time 
lattice A, i.e. (x, t) G A. This way, we impose a UV cutoff that eliminates the details of 
those processes occurring deep in the dissipative regime. Then, the measure is given by 
[du] — ^ n{a; t) e a'^^(^' ^^'^ action in (E]) needs to be discretized appropriately. We 
replace the dynamics ([1]) with a finite-difference equation with backward-time discretization 



dfU + udxU —J- -{u{t) — u(t — e)) -|- u{t — e) d^u^t — e) 



(6) 



where e is the lattice spacing in time direction. This ensures the correct dynamics in the 
continuum limit 23| . For the advective term we take the anti-symmetric spatial derivative 



dxU — 7- —{u{x + a) — u{x — a)) 



(7) 



where a is the lattice spacing in the spatial direction. With this choice of discretization 
the problem is amenable to a local over- relaxation algorithm 2J|. Starting from an initial 
configuration {u{x , t) , {x , t) G A} the set of single-site variables is updated iteratively by 



the successive application of a transition probability P(u(x, t) — > u'{x,t)). We use the high- 
quality ranlux (pseudo) random number generator 25| which is essential for large-scale 
lattice simulations. Specific improvements, e.g. Chebyshev acceleration j26| significantly 
reduce thermalization and autocorrelation times for the relevant observables in the inertial 
range. 

It is necessary to map the discretized theory to its continuum counterpart and one has 
to ensure that the parameters are well-defined in the continuum limit. For that purpose the 
kinematic viscosity is identified with u = vo? /e where v is the viscosity in lattice units, and 
the Reynolds number scales as Re oc . Furthermore, we have to ensure that the relevant 
scales of the system are resolved. In particular, we have to ensure that the dissipation scale 
fits on the lattice, i.e. rj = Re~^/^L > a where L is the IR scale present in our system as a 
consequence of the finite lattice size. One may immediately recognize that this imposes a 
hard constraint on the realization of lattice simulations - fully developed turbulence requires 
a large computational effort where the number of lattice sites in the spatial direction scales 
as oc Re^/^ for given L. In practice, we are therefore bound to work at non-zero viscosity u. 

Simulations at moderate to high Reynolds numbers require massively parallel architec- 
tures. State of the art simulations at Re = 64 and lattice size 254 x 1024 (space x time 
direction) run on up to 512 processors in parallel. Structure functions are evaluated over an 
ensemble of configurations generated by the MCMC algorithm as described in the previous 
paragraphs. For every configuration we measure velocity differences from a randomly cho- 
sen starting point. This dramatically reduces autocorrelations for our sample. In Fig. [2^ 
we show, as an example, the n = 5 order structure function calculated for an ensemble 
from a 254 x 1024 lattice simulation, averaged over nearly 5 x 10^ statistically independent 



field configurations. To determine t 
problem in the literature (see e.g. 



le scaling range a priori is difficult, and a well-known 
l|). Here, we employ a working definition where it is 
defined as the range of scales that minimizes the of ^ linear least-squares (LLS) fit to the 
fifth order structure function in the log-log plot. The corresponding region is indicated in 
Fig. [2ti. For comparison we have included the values of the local slope (evaluated over three 
consecutive points) in the inset. We identify a plateau where the local exponents are nearly 
constant - this defines the value of the scaling exponent. In general, with this method, we 
cannot rule out subleading terms or possible logarithmic corrections that may influence the 



scaling behavior 
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28| . In fact, such a situation is very likely and can lead to the 
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FIG. 2: (a) Log-log plot of the structure function of order n = 5 with a linear scaling function 
plotted for comparison. Vertical bars indicate the region for the extraction of scaling exponents. 
Inset shows the local slopes versus r. (b) Structure function scaling exponents Cn versus order n. 
The black curve indicates a bifractal scaling behavior. 



appearance of multiscaling 18|. While in principle these contributions should be taken into 
account for the accurate determination of the scaling behavior, in practice it is difficult to 
distinguish different types of scaling contributions without any further assumptions. We 
obtain the scaling spectrum (Fig. |2b) where the error bars given are those of the LLS fit in 
the scaling range. Clearly, the n = 5 data point in Fig. |2]d has minimal error which follows 
simply from our definition of the scaling range. We see that the scaling exponents are close 
io the bifractal scaling prediction js], [oj, and within error bars agrees with the results of 
18|, obtained at high spectral resolution. As a last remark we want to add that we have 
not used extended self-similarity (ESS) {29 1 at any point in our analysis. Though ESS may 
enlarge the effective scaling range we found that it can suggest a clean scaling behavior even 
if sub leading terms are present. 

Since we are dealing with a finite system both in space and time one may expect boundary 
effects. In our simulations we have chosen periodic boundary conditions in space and fixed 
(Dirichlet) boundary conditions in time. This way we eliminate zero mode effects from the 
dynamics. For a space-time lattice of infinite extent in the time direction the probability 
measure ([5]) defines a stationary process, i.e. correlation functions will only depend on 
time differences. We have checked this property explicitly in our analysis - far from the 
boundaries, in the bulk of the configurations, the system is approximately in a stationary 
state. 
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FIG. 3: Probability distribution functions 7^(An, r) as a function of the dimensionless variable 
(f) = An/[(Ati^)]^/^ plotted for different values of r. (a) Collapse of the PDF in the universal 
regime (blue). In the energy-containing range (red) the fluctuations become Gaussian - the random 
forcing dominates ~ whereas in the dissipative regime (orange) fluctuations are strongly enhanced, 
(b) Scaling region for the left tail of the PDF. The black line indicates the scaling prediction with 
exponent 7 = — 4. 

In the continuum both the action in (|5]) and the measure are invariant under the set of 
Galilean transformations 



X 



X + r , u{x, t) — )■ u{x + r,t) + V , r = vt . 



(8) 



To avoid an over counting of physically equivalent field configurations one should eliminate 
these modes by the Faddeev- Popov procedure (23[. While gauge fixing is unavoidable for 



generic correlators 
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31| this is not so for velocity differences, as solely considered in this 
work which are manifestly invariant under transformations ([8]). 

One may also check the statistics for velocity differences directly on the level of the 
probability distribution functions V{Au,r). This gives valuable qualitative information on 
the physical behavior in our simulations of Burgers turbulence. In Fig. |3t we show the PDF 
of velocity differences for a set of values of the separation r, where we use the dimensionless 
variable = Am/[(Am^)]^/^ to quantify the fluctuations. At large scales, far from the inertial 
range we clearly recognize the effects of the random Gaussian forcing (red). In the dissipative 
region the left tail of the PDF is especially pronounced and captures the strong fluctuations 
described by the shocks (orange). For separations 77 ^ r ^ L in the inertial range we see 
that the PDF V{Au,r), plotted for three different values of r, nicely collapse onto each 
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other (blue). In particular, in the regime where the fluctuations are much smaller than 
the root-mean-square velocity | An| ^ Urms, the PDF of velocity differences has a universal 
scaling form 

V{Au,r) =r-'f{Au/r') , (9) 



where z is the dynamic exponent [32|. In the asymptotic region —Au/r^ ^ 1 where Am < 
we expect the algebraic scaling V{Au,r) oc (Au)'^ with exponent 7 = —4. The relevant 
region is shown in Fig. (indicated by the arrow) and Fig. 13b- The corresponding scaling 
prediction with exponent 7 = —4 is plotted for comparison as the black line in Fig. [3b. 
Though our statistics are not sufficient to give a tight prediction on the scaling exponent, 
indications for the conjectured scaling behavior can be inferred from Fig. |3]d. 

At this point we want to give a short remark on some issues that arise when turning 
to incompressible three-dimensional Navier-Stokes turbulence. It is well-known, that the 
inclusion of the pressure term is one of the main obstacles in simulations of turbulence, 
as the requirement of incompressibility introduces non-local interactions. In the functional 
integral formulation this leads to a non-vanishing Fadeev-Popov determinant that can be 



treated by standard procedures (see e.g. |33|). In our lattice approach, we have chosen to use 
a local update over-relaxation algorithm that proved to be quite efficient for our purposes. 
The long-range correlations in our system imposed by the forcing however, prohibit any 
attempt to parallelize in the spatial direction. This poses a severe problem when turning 
to higher dimensions. For that purpose it is absolutely mandatory to switch to a global 



update algorithm as, e.g. a Hybrid Monte Carlo algorithm (3J] that is usually used in 
standard lattice calculations when non-local actions are considered. The implementation of 
non-trivial spatial boundary conditions is an interesting possibility when turning to higher 
dimensions. In principle, there are no restrictions on the type of boundary conditions in our 
simulations. However, in such a case the two-point correlation function of the forcing will 
not be restricted by symmetry arguments and may take a rather difficult form. 

The Burgers equation provides an ideal setting to understand systematic effects at finite 
Reynolds numbers and finite lattice sizes in the framework of Monte Carlo simulations. We 
have demonstrated that our simulations are able to reproduce the well-known anomalous 
scaling behavior in the Burgers model. It is important to remark that this is possible without 
exploiting the integrability property of the Burgers equation, as was done e.g. with a fast 



Legendre transform algorithm in 
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35| . Thus, Monte Carlo simulations are directly appli 



cable to other physical systems of interest where it is important to have alternative methods 
available to determine the statistical behavior of small-scale fluctuations. Certainly, the 
numerical efficiency of our method is an issue, and there is room for improvement. Specifi- 
cally, one may ask if a global Hybrid Monte Carlo algorithm will improve the performance 
of our simulations. Another question isolated from the numerical efficiency relates to the 
evaluation of the scaling spectrum and the role of subleading scaling corrections. For that 
matter it is essential to distinguish possibly different types of scaling contributions. Further 
work in this direction is in progress. This work presents a first important step towards 
the determination of scaling behavior for the case of full three-dimensional Navier-Stokes 
turbulence in the framework of Monte Carlo simulations. 
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